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Abstract 

Background: The major histocompatibility complex (MHC) genes are one of the most important genetic systems in 
the vertebrate immune response. The diversity of MHC genes may directly influence the survival of individuals 
against infectious disease. However, there has been no investigation of MHC diversity in the Asiatic black bear 
{Ursus thibetanus). Here, we analyzed 270-bp nucleotide sequences of the entire exon 2 region of the MHC DQB 
gene by using 188 samples from the Japanese black bear {Ursus thibetanus japonicus) from 12 local populations. 

Results: Among 185 of 188 samples, we identified 44 MHC variants that encoded 31 different amino acid 
sequences (allotypes) and one putative pseudogene. The phylogenetic analysis suggests that MHC variants 
detected from the Japanese black bear are derived from the DQB locus. One of the 31 DQB allotypes, Urth-DQB*01, 
was found to be common to all local populations. Moreover, this allotype was shared between the black bear on 
the Asian continent and the Japanese black bear, suggesting that Urth-DQB*01 might have been maintained in the 
ancestral black bear population for at least 300,000 years. Our findings, from calculating the ratio of 
non-synonymous to synonymous substitutions, indicate that balancing selection has maintained genetic variation 
of peptide-binding residues at the DQB locus of the Japanese black bear. From examination of genotype frequencies 
among local populations, we observed a considerably lower level of observed heterozygosity than expected. 

Conclusions: The low level of observed heterozygosity suggests that genetic drift reduced DQB diversity in the 
Japanese black bear due to a bottleneck event at the population or species level. The decline of DQB diversity 
might have been accelerated by the loss of rare variants that have been maintained by negative 
frequency-dependent selection. Nevertheless, DQB diversity of the black bear appears to be relatively high 
compared with some other endangered mammalian species. This result suggests that the Japanese black bears 
may also retain more potential resistance against pathogens than other endangered mammalian species. To 
prevent further decline of potential resistance against pathogens, a conservation policy for the Japanese black 
bear should be designed to maintain MHC rare variants in each local population. 

Keywords: Balancing selection, Conservation genetics, DQB, Genetic diversity, Japanese black bear, Major 
histocompatibility complex, Ursidae, Ursus thibetanus 



Background 

In Japan, the Asiatic black bear (Ursus thibetanus) is 
classified as the subspecies {Ursus thibetanus japonicus). 
The Japanese black bear currently inhabits two main 
islands in Japan, Honshu and Shikoku. Historically, these 
bears were also distributed in Kyushu, but no native 
bears have been captured there for several decades [1,2]. 
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However, there have been several reports of sightings of 
an animal believed to be a black bear in Kyushu [3]. Re- 
search on the genetic diversity in the black bear has 
begun in recent years and has been developed with the 
use of genetic markers such as mitochondrial DNA 
(mtDNA) sequences and microsatellites [4-12]. Using 
mtDNA markers we have found that the Japanese black 
bear is genetically quite distinct from other Asiatic black 
bears [8]. Sequences of the mtDNA control region indi- 
cated that the Japanese black bear could be divided gen- 
etically into three groups: western, eastern, and Shikoku/ 
Kii-hanto populations [7,8]. Analyses of mtDNA and 
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microsatellite markers [6,8] revealed lower genetic diver- 
sity and higher genetic differentiation in the western 
population than in the other populations. In fact, in 
most of the management units in western Japan the spe- 
cies has been designated as "endangered" in the Red 
Data Book of Japan [13] because of the fragmentation 
and isolation of their habitats. 

The major histocompatibility complex (MHC) genes 
play an important role in the resistance to various 
pathogens, especially in wild animals. If MHC poly- 
morphism declines, resistance to infectious diseases will 
decrease [14]. The MHC genes are thought to influence 
the survival of species or populations against pathogens. 
Therefore, MHC genetic variability in local populations 
is regarded as a potentially important index in conserva- 
tion genetics. 

In the human genome, the MHC region is located on 
the short arm of chromosome 6 at 6p21.3 and covers 
roughly 4 Mbp comprising 224 genes [15]. The human 
MHC region encodes at least three functionally and 
structurally different types of molecules, class I, class II, 
and class III. The MHC class II region is composed of 
five sub-regions (DP, DQ, DR, DM, and DO). Each sub- 
region contains at least one set of A and B genes, and 
these encode two polypeptide chains, a and chains, re- 
spectively, which make up a heterodimeric molecule. For 
example, in a DQ molecule, the DQA gene encodes the 
a-chain, while DQB encodes the /?-chain. The DP, DQ, 
and DR molecules are primarily important for the pres- 
entation of peptides and they are classified as classical 
class II molecules. 

For the last three decades, studies of MHC allele varia- 
tions have been performed not only in humans and la- 
boratory animals but also in populations of wild 
mammals. Many such studies of wild mammals have 
demonstrated a relationship between low levels of MHC 
diversity and decreases in species population fitness (see 
below). MHC diversity was examined in pinniped spe- 
cies, which are believed to be closely related phylogenet- 
ically to ursine species [16,17]. Hoelzel and colleagues 
[18] reported that genetic variation at the DQB locus in 
the southern elephant seal {Mirounga leonina) did not 
appear to be low, whereas the northern elephant seal 
{Mirounga angustirostris), which was thought to have 
experienced a severe population bottleneck, exhibited 
much less variation at this locus. Lento et al. [19] and 
Weber et al. [20] demonstrated low genetic diversities at 
the MHC loci in the New Zealand sea lion (Phocarctos 
hooker i), and the northern elephant seal. The census 
sizes of these pinniped species have also diminished in 
the recent past probably due to human activity. In the 
case of non-carnivore mammals, heterozygosity of the 
DRB gene in a population of the striped mouse (Rhabd- 
omys pumilio), has been shown to influence infection 



status [21]. A low level of diversity at the DAB locus (an 
MHC class II gene in marsupials) in the western barred 
bandicoot (Perameles bougainville) may influence its vul- 
nerability against novel pathogens [22]. 

Many studies have revealed a correlation between 
MHC diversity and resistance against infectious disease; 
however, others have demonstrated that this is not the 
case universally. For example, genetic variation in DQ 
and DR sequences in California sea lions (Zalophus cali- 
fornianus), which is a thriving species, is lower than that 
in several other mammalian species [23,24]. The Tas- 
manian devil (Sarcophilus harrisii), which is an endan- 
gered species, has low MHC diversity, but individuals 
with restricted MHC repertoires (These individuals pos- 
sessing only one of two phylogenetic groups) may be re- 
sistant to the devil facial tumor disease [25]. However, 
their MHC variants were not assigned to loci because of 
a lack of genomic information. Thereafter, the research 
group has shown that there is no direct evidence that 
MHC genetic difference between tumor-free and 
infected animals is associated with the resistance against 
the infection although certain genotypes at two MHC- 
linked microsatellite markers seem to be more fre- 
quently detected from healthy animals [26]. The DRB 
locus in the mountain goat (Oreamnos americanus) 
exhibited low genetic diversity, but this does not appear 
to have affected the resistance to infectious diseases in 
this species [27]. 

For Ursidae, studies of MHC diversity are somewhat 
limited. In the giant panda {Ailuropoda melanoleuca), 
MHC studies are relatively advanced and have been used 
to propose some possible implications for captive man- 
agement [28-34]. Low genetic diversity at the Aime-DRB 
and Aime-DQA loci has been observed in the giant 
panda [28,30]. Wan et al. [31] attempted to understand 
the genomic structure and evolutionary history of Aime- 
MHC genes. The study has showed that the giant panda 
possesses a single DYB-like pseudogene (the DY sub- 
region was thought to be ruminant-specific MHC genes, 
but may be shared by members of Laurasiatheria super- 
order.). Chen et al. [34] examined genetic diversity at 
seven Aime-MHC class II genes from 121 individuals 
and showed the presence of two monomorphic and five 
polymorphic loci. In brown bears (Ursus arctos) three 
DQA variants and 19 DRB variants were identified in 32 
and 38 bears, respectively, although loci have not been 
assigned [35,36]. Kuduk and colleagues reported the ex- 
pression of three MHC class I, two DRB, two DQB and 
one DQA loci in the brown bear [37]. Other than these 
studies on the giant panda and brown bear, there have 
been no published reports on the MHC genetic diversity 
in the Ursidae. 

Although the direct contribution of MHC genes to the 
immunological fitness of animals is still unclear, the 
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relatively high genetic diversity of the MHC reflects the 
abundance in repertoire of binding peptides, derived 
from pathogens. We thus believe that MHC diversity po- 
tentially contributes to resistance against environmental 
pathogens and species survival, indicating that study of 
these loci is essential for designing a plan for conserva- 
tion management of animals. Because of large variability 
of MHC, the examination of MHC diversity is also use- 
ful for the management of endangered wild-animal 
populations that show no variation in genetic markers 
such as microsatellite loci [38]. In our previous study, 
we reported a functional DQB sequence (1,150 bp) in 
the Japanese black bear, identified using the rapid ampli- 
fication of cDNA ends method [39]. Here, we investi- 
gated polymorphism in exon 2 of DQB that contains the 
hypervariable region involved in binding antigens for 
the Japanese black bear to consider the potential impli- 
cations for resistance against pathogens and estimate the 



status of local populations on the basis of their MHC 
variation. 

Results 

Features of the DQB exon 2 variants in the Japanese 
black bear 

The entire exon 2 region of the DQB locus (270 bp) was 
determined for the Japanese black bear. The amplified 
region corresponded to amino acids 6 to 94 of the /? 
chain of the human DR molecule [40,41]. The frequency 
of variant typing success was lowest for tooth samples 
(Ca. 54%), and was highest for blood (Ca. 92%). This re- 
sult probably reflects the difference in the amount of 
gDNA in each sample type. From 188 bear samples in 
12 of 19 conservation and management units designed 
by Yoneda [42] (Figure 1), we detected 47 DQB variants. 
However, the application of criteria for genotyping used 
in this study (see Methods) reduced the number to 44 
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Figure 1 Geographical distribution of DQB allotypes and haplotypes of mtDNA control region. The inner pie chart indicates the observed 
frequencies of DQB allotypes in the Japanese black bear. The outer pie chart indicates and the observed haplotype frequencies of the mtDNA control 
region [8]. Twelve local populations based on the conservation and management units for the Japanese black bear [42] are shown as bold circles. The 
number of samples for DQB allotypes from each location is depicted in the center of each corresponding chart, and the number in parenthesis 
indicates the sample size used for mtDNA analysis of our previous study [8]. Superscript a represents the number of samples analyzed by Ishibashi and 
Saitoh [5] in addition to our previous study [8]. Superscript b indicates the number was obtained by identification of individuals using microsatellite 
analysis [90]. Superscript c indicates the putative pseudogene. Superscript d indicates haplotypes of mtDNA identified by Ishibashi and Saitoh [5]. 
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variants [DDBJ: AB634514-AB634536] and to 185 bear 
samples (Figures 1, 2 and Table 1). In the present study, 
no more than two different variants per individual were 
detected from the clones of the PCR products. 

Thirty-one distinct amino acid sequences (Urth- 
DQB*01 to Urth-DQB*31) were translated from nucleo- 
tide sequences of the 43 variants (Urth-DQB*0101 to 
Urth-DQB*3101) (Figure 3). Because one variant exam- 
ined (Urth-DQB*3201) included a premature stop codon 
in the sequence, it was recognized as a pseudogene 
(Figure 3A). This Urth-DQB*3201 was detected from 



one individual which was heterozygous with this variant 
and Urth-DQB*0101. An NCBI BLAST search showed 
that sequence similarities between bear Urth-DQB and 
panda Aime-DQB variants were 84% to 91% and 93% to 
96% similarity at the amino acid and nucleotide levels, 
respectively [34]. Of all the variants, Urth~DQB*0101, 
*0104, *0401, *1001, and *2601 correspond to the Urth- 
DQB*a to Urth-DQB*e variants detected in our previous 
study [39], at the DNA level. Thus, in the present study 
we additionally identified 39 novel variants at the DNA 
level in 178 samples. The nucleotide sequence of Urth- 
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Figure 2 Alignment of segregating sites among the DQB exon 2 variants in the Japanese black bear. Dots indicate identity with the 
nucleotides of Urth-DQB*0101 . The first two numeric characters after an asterisk in a variant name represent an allotype. 
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Table 1 The observed number of nucleotide DQB variants detected from each conservation and management unit 

Conservation and CS GAI SO EM SA NCA HO NK EC WC KH SK Total 
management unit 

No. of sample 5 21 1 13 1 12 36 23 8 53 3 

No. of variant 2 4 1 7 1 10 18 8 7 10 3 
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Table 1 The observed number of nucleotide DQB variants detected from each conservation and management unit 

(Continued) 
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Urth-DQB*3101 
Urth-DQB*3201 



CS: Chokai Sanchi; GAI: Gassan-Asahi lide; SO: Southern Ouu; EM: Echigo-Mikuni; SA: Southern Alps; NCA: North-Central Alps; HO: Hakusan-Okumino; NK: Northern 
Kinki; EC: Eastern Chugoku; WC: Western Chugoku; KH: Kii Hanto; SK: Shikoku. 



DQB*0101 was identical to that found in exon 2 of a 
cDNA sequence derived from mRNA of a bear sample 
in our previous study [39]. This suggests the amplifica- 
tion of a functional DQB locus in the present study, al- 
though there is no evidence of protein translation. 
Deletions and insertions were not detected among the 
sequences. There were 58 segregating sites among the 44 
different nucleotide sequences (Figure 2). Amino acid 
substitutions were observed at 30 sites among the 31 



allotypes (Figure 3). While 11 of these substituted sites 
were located in the putative peptide-binding residues 
(PBRs) [40], 19 corresponded to non-PBR positions. The 
presence of more substitution sites in the non-PBRs 
appears to suggest a loss of the open reading frame in the 
bear DQB gene. However, there were no frameshift or 
nonsense mutations in Urth-DQB sequences with the ex- 
ception of the Urth-DQB ^3201 although our study was 
limited to exon 2 (Figure 3). In addition, amino acids at 8 
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1 I Positions of PBRs estimated based on the human DR molecule (Brown et a/.1993) 
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Figure 3 Amino acid sequences and the variability in the exon 2 of bear DQBs. (A - B) Gray boxes indicate positions of the putative 
peptide binding residues (PBRs) deduced from the amino acid position in the human DR molecule, and the numbers across the sequence 
correspond to amino acid positions based on the chain of the human DR structure [40]. A filled asterisk represents the putative site under 
balancing selection identified by both of PAML and OmegaMap programs, while an open asterisk is that by the OmegaMap only. (A) Alignment 
of the predicted amino acid sequence translated from nucleotides of MHC DOB variants in the Japanese black bear. Dots indicate identity with 
amino acids of Urth-DQB*01. Urth-DQB*32 includes the putative premature stop codon within the sequence. (B) The variability level at amino 
acid residues among DQB variants of the Japanese black bear. The variability level was estimated by a Wu-Kabat plot [85]. 
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of 12 PBR codons that showed a monomorphic pattern 
were also conserved in the human DQB1 locus, and most 
of the segregating sites at the non-PBRs among bear DQB 
variants were singletons. Although we cannot completely 
exclude the possibility that pseudogenes are present in 
our data set, the fact that our results show more substitu- 
tion sites in the non-PBRs than PBRs indicates that it is 



unlikely that we amplified pseudogenes. Therefore, we 
defined all bear DQB variants as the putative functional 
DQB molecule with the exception of Urth-DQB*3201. 
Thirty-one DQB allotypes could be classified into 15 
types of PBR (i.e., 15 different amino acid sequences of 
the PBRs; Figure 4). In 184 bears (excluding the individual 
with the putative pseudogene Urth~DQB*3201), 57 of 64 



(A) 



70 



Urth-DQB *31 

Urth-DQB *25 

Urth-DQB *23 
Urth-DQB* 

Urth-DQB* 15 

Urth-DQB* 10 

Urth-DQB *05 

Urth-DQB *01 
Urth-DQB *1 



- Urth-DQB *29 
Urth-DQB *24 

■ Urth-DQB* 12 

77]_l Urth-DQB *21 
701 Urth-DQB *22 
Urth-DQB *30 

Urth-DQB *03 
Urth-DQB* 17 
Urth-DQB* 19 
Urth-DQB *28 
Urth-DQB *09 
Urth-DQB* 14 



Urth-DQB *07 
Urth-DQB* 11 
Urth-DQB *20 

Urth-DQB *02 

Urth-DQB *04 
Urth-DQB *06 
Urth-DQB *08 
60j Urth-DQB *13 

Urth-DQB *26 



Urth-DQB *27 



- Giant pa 

— Giant r 



da (Aime-DQB1*01: GQ496182) 
nda(Aime-DQB2*01: GQ496188) 




Aime-DQB1*01 (GQ496188) 
(GQ496182) 

Figure 4 Neighbor-joining tree and network based on amino acids in the PBR of DQB allotypes. (A) A neighbor-joining (NJ) tree was 
constructed by using the number of differences method. DQB sequences of the giant panda are used as outgroups. Only bootstrap values over 
50% are shown in this figure. (B) Neighbor-Net network is also based on the sequences of NJ tree. 

V . ) 
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heterozygotes possessed a combination of different types 
of PBR sequences. 

On the basis of the nucleotide sequences of the 
black bear DQB variants, we performed a search for 
intragenic recombination by using the GENECONV 
program [43]. Permutation testing indicated the possi- 
bility of recombination between some pairs of variants 
(Urth-DQB*0101, Urth-DQB*0104, Urth-DQB*0109, Urth- 
DQBW402, Urth-DQB*0601, Urth-DQB*0701, Urth-DQB* 
0901, Urtii-DQB*1101, Urth-DQB* 1301, Urth-DQB*2001, 
Urth-DQB*2S01, Urth-DQB*2701, Urth-DQB*2901, Urth- 
DQB*3001, Urth-DQB*3101, and Urth-DQB*3201). How- 
ever, the Bonferroni-corrected BLAST-like test in 
GENECONV indicated no significant values (P = 1.0) 
for recombination events. In addition, not only the 
method described by Satta [44] but also the RDP4 pro- 
gram [45] indicated no significant signals for recombin- 
ation or gene conversion. Thus, we used all the DQB 
variants of the Japanese black bear for further analyses. 

Gene tree based on exon 2 sequences of ursine MHC 
class II genes 

A maximum likelihood (ML) tree based on amino acid 
sequences of partial DQB and DRB exon 2 from the 
Ursidae species was constructed (Figure 5). In this ML 
tree, the DQB sequences formed a clearly separated clus- 
ter from the DRB sequences, in agreement with Yasuko- 
chi et al. [39]. These trees showed that the 31 sequences 
from the Japanese black bear formed a monophyletic 
group with the DQBs of the giant panda. The Bayesian 
and Neighbor- joining (NJ) [46] trees also showed a 
monophyletic clade of the black bear DQBs and panda 
DQBs (see Additional File 1: Figure SI). 

In the ML tree, DRBs of some ursine species were 
mixed with each other although the orthology of their 
sequences was unclear, indicating a possible trans- 
specific mode in the DRBs of ursine species. In the 
black bear/panda DQB clade, twenty black bear DQB 
allotypes (Urth-DQB*02 to *04, *06 to *09, *11 to *14, 
•17, *19 to *22, *26 to *28, and *30) appear to cluster 
more closely with panda rather than other black bear 
DQBs, suggesting the trans-species polymorphism. 
However, the bootstrap value is low (41%), and the 
pairwise distance among the twenty bears and other 
bear DQB variants showed closer relationships than 
those between black bear and panda DQBs. For ex- 
ample, the pairwise distances between the Urth- 
DQB*04 and other black bear DQBs ranged from 0.01 
to 0.08 (mean value, 0.06), while those between the 
Urth-DQB*04 and panda DQBs ranged from 0.09 to 
0.16 (mean value, 0.13). In addition, the Bayesian tree 
showed that DQB sequences of Japanese black bear 
formed a distinct clade from those of giant panda (see 
Additional File 1: Figure SI). Therefore, the allelic 



lineage of twenty black bear DQBs and panda DQBs 
may be different even though the ML tree appears to 
show the pattern of trans-species evolution. 

Selection of the DQB exon 2 region 

We estimated the mean numbers of synonymous and 
nonsynonymous substitutions within the 267-bp re- 
gion of DQB exon 2 variants in the Japanese black 
bear by using a modified Nei-Gojobori method [47] 
(Table 2). For the black bear, the ratio of nonsynon- 
ymous to synonymous substitutions (d N /d s ) was ap- 
proximately 3.83 in the PBRs, 0.63 in the non-PBRs, 
and 1.48 in the entire exon 2 region. Although three 
statistical tests did not indicate recombination among 
all the DQB variants, according to the BLAST-like test 
without the Bonferroni-correction in GENECONV 
program (see above), 16 variants were shown to have 
possibly undergone recombination or gene conversion. 
Recombination would change the d N /d s value, how- 
ever, even if these possible recombinants were 
removed, ratios of d N /d s showed slightly higher 
values: 4.20 in the PBRs, 0.58 in the non-PBRs, and 
1.52 in the entire exon 2 region. This shows that in- 
clusion of the possible recombinants does not affect 
the result. To evaluate the d N /d s of the black bear, we 
compared them with those of other mammalian DQB 
or DRB loci (Table 2). In the family Ursidae, the d N /d s 
ratio within the PBRs of DQB in the black bear was 
higher than that of the giant panda DQBs, while it was 
lower than that of brown bear DRBs (but these DRB 
variants are derived from at least two loci). Compared 
with other mammals, the d N /d s value of the black bear 
did not appear to be low. In addition, by using the 
PAML program with models that describe various se- 
lection pressures (M0 to M8) [48-50], ratios of co were 
estimated for black bears (Table 3). Values of co in the 
black bear were significantly higher than unity (log like- 
lihood ratio test) under M2a and M8 models that allow 
for balancing selection (co 2 = 7.0 in M2a and co = 5.7 in 
M8). The likelihood ratio test of hypotheses (M2a vs 
Mia and M8 vs M7) showed statistically significant 
values (P < 0.001), suggesting a balancing selection model 
is a better fit to the data compared to other models. Thus, 
the present results indicate that exon 2 of DQB in the black 
bear has been a target of balancing selection. The Bayesian 
analysis also showed that five codons were potentially 
under balancing selection. The mean value of co, estimated 
with Bayesian inference by using the OmegaMap program 
[51], was also larger than unity (co = 1.8). This analysis also 
suggested that 11 codons were probably under balancing 
selection. Of these 11 codons, five were estimated by both 
the PAML and the OmegaMap programs. A majority of 
the selected sites (|39, |328, (330, and (337) were located in 
the predicted PBRs (Figure 3). 
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?Giant panda Aime-DRB3*05 (GQ496169) 
- Giant panda Aime-DRB 3*06 (GQ496170) 
Giant panda Aime-DRB 3*03 (GQ496167) 
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Polar bear Urma-DRB*22 (AF458935) 
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Polar bear Urma-DRB*20 (AF458933) 
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Polar beai- Urma-DRB*08 (AF458921) 
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)J — Polar bear Urma-DRB*16 (AF458929) 



Polar bear Urma-DRB*25 (AF458938 
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- Brown bear Urar-DRB*03 (AB490459) 
-Sun bearUrml-DRB*01 (AB490483) 
-American black bear Uram-DRB*01 (AB490480) 

Sun bear Urml-DRB*02 (AB490484) 

-Brown bear Urar-DRB *01 (AB490457) 

American black bear Uram-DRB*02 (AB490481) 
- American black bear Uram-DRB*03 (AB490482) 
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in Sloth bear Urur-DRB*02 (AB490486) 
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r Giant panda Aime-DRB *5 (AY895 159) 
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m Giant panda Aime-DRB *7 (AY895161) 
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Figure 5 Maximum likelihood tree for the exon 2 region of MHC class II genes. This tree is constructed based on 79 amino acids in the 
MHC class II DRB and DOB exon 2 region. The genetic distances are computed by the model with gamma distribution. Only bootstrap values 
over 50% are shown in this figure. Numbers in parentheses are Genbank accession numbers. 
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Table 2 Comparison of non-synonymous (c/ N ) and synonymous (c/ s ) substitution rates at the DQB exon 2 region 



Species 


Gene 


No. of 
sample 


No. of 
variant 


Length 


PBRs 








Non-PBRs 






Entire exon 2 region 






d N /d s 




d u 




d H ld s 








d H ld s 




d H 




d s 




Ursidae 






































Japanese black bear 3 


DOB 


184^ 


43* 


267 bp 


3.83 


0.11 


± 


0.04 


0.63 


0.02 


± 


0.01 


1.48 


0.04 


± 


0.01 


0.03 ± 


0.01 


Giant panda b 


DQB1 


121 


6 


267 bp 


1.33 


0.02 


± 


0.06 


0.47 


0.01 


± 


0.02 


0.71 


0.02 


± 


0.05 


0.01 ± 


0.03 


Brown bear* c 


DRB 


38 


19 


267 bp 


4.09 


0.06 


± 


0.29 


0.72 


0.02 


± 


0.05 


1.53 


0.02 


± 


0.07 


0.02 ± 


0.11 


Other carnivores 






































African wild dog d 


DRB 


320 


19 


267 bp 


5.28 


0.06 


± 


0.21 


0.56 


0.01 


± 


0.02 


1.72 


0.02 


± 


0.04 


0.02 ± 


0.07 


Namibian cheetah* 6 


DRB 


139 


4 


237 bp 


2.01 


0.06 


± 


0.25 


0.37 


0.01 


± 


0.05 


0.76 


0.04 


± 


0.13 


0.02 ± 


0.10 


European mink f 


DRB 


20 


9 


228 bp 


1.78 


0.06 


± 


0.16 


0.26 


0.01 


± 


0.02 


0.65 


0.03 


± 


0.09 


0.02 ± 


0.06 


Non-carnivores 






































Water buffalo* 9 ' h 


DOB 


210 


22 


243 bp 


2.65 


0.07 


± 


0.37 


0.92 


0.03 


± 


0.13 


1.34 


0.03 


± 


0.14 


0.03 ± 


0.19 


Ma's night monkey 1 


D0B1 


19 


14 


246 bp 


2.71 


0.05 


± 


0.19 


0.95 


0.02 


± 


0.07 


1.39 


0.02 


± 


0.07 


0.02 ± 


0.10 


Yellow necked mouse j 


DRB 


146 


27 


216 bp 


7.75 


0.09 


± 


0.38 


1.13 


0.02 


± 


0.06 


2.46 


0.02 


± 


0.05 


0.02 ± 


0.12 



The c/ N and d s are estimated by the modified Nei-Gojobori method. 

a : present study; b : Chen et al. [34]; c : Goda et al. [36]; d : Marsden et al. [91]; 

e : Castro-Prieto et al. [92]; f : Becker et al. [93]; 9 : Niranjan et al. [94]; h : Sena et al. [95]; 

': Diaz et al. [96]; j : Meyer-Lucht and Sommer [97]; k : An individual that has the putative pseudogene, Urth-DQB*3201 , is excluded. 
*: The number of MHC variants may be derived from multiple loci. 



Estimation of the divergence time of DQB variants in the 
Japanese black bear 

To estimate the divergence time of DQB variants in the 
Japanese black bear, we calculated the maximum genetic 
distance at synonymous sites that corresponded to the 
time to most recent common ancestor (TMRCA). The 
resulting value was 0.08 ± 0.03, and, by assuming 3.5 x 
10~ 9 per site per year as the average synonymous substi- 
tution rate for mammalian nuclear DNA [52], the diver- 
gence time of the bear DQB variants was estimated to be 
about 11 million years ago (MYA). Even if the possible 
recombinants described above were removed, the result 
did not change. 



Geographical distribution and genetic variation in the 
MHC exon 2 variants in the conservation and 
management units 

Five DQB allotypes (Urth-DQB*01, Urth-DQB*02, 
Urth-DQB*03, Urth-DQB*04, and Urth-DQB*07) were 
shared among the management units (Figure 1 and 
Table 1). In particular, Urth-DQB*01 was observed 
throughout all the management units. This allotype, 
Urth-DQB*01, corresponds to six nucleotide DQB var- 
iants (Urth-DQB*0101 to Urth-DQB*0109). The variant 
Urth-DQB*0101 was detected in a particularly large 
number (133) of bears. This suggests that Urth-DQB*01 
may have been dominantly maintained in the black bear 



Table 3 Inference of balancing selection for 43 black bear DQB variants 

Method Model code Parameter estimates Likelihood or confidential interval Positively selected sites 

Maximum likelihood (PAML) M0 (one ratio) go = 0.41 -1570 



Mia (nearly neutral) 


p 0 = 0.85,0! = 0.15 


-1386 








M2a (balancing selection) 


co 2 = 7.00 


-1326 


P9** 


, (326**, P28**, 
P37** 


P30**, 


M3 (discrete) 


co 2 = 7.98 


-1325 








M7 (beta) 


p = 0.05, q = 0.07 


-1394 








M8 (beta and go) 


co = 5.71 


-1337 


P9** 


(326**, P28**, 
P37** 


P30**, 



Bayesian inference (OmegaMap) - a) = 1 .80 a u: 0.30 - 1 2.49 a (95% HPD) P9** 7 (322**, (326**, P28*, 

(329**, P30**, P37*, (350*, 
(357*, P68*, (384* 

The putative pseudogene (Urth-DQB*3201) is excluded from this analysis, co = d N I d s . p 0 : the proportion of conserved sites with 0 < oj < 1. (=1 - p 0 ): the 
proportion of neutral sites. Values of p and q are parameters of the beta distribution. *: posterior probability, P < 0.95; **: posterior probability, P < 0.99. a : a mean 
value. Bold characters indicate positions of the putative peptide binding residues (PBRs) deduced from the amino acid position in the human DR molecule [40]. 
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population in recent times. This is probably because of 
its important immune role against pathogens or for 
demographic reasons. 

On the basis of the variant frequency in each of the six 
management units from which more than 10 individuals 
were sampled, the expected heterozygosity (H e ) ranged 
from 0.33 to 0.85 (mean H e , 0.65) (Table 4). Values of 
observed heterozygosity (H Q ) ranged from 0.19 to 0.69 
(mean H Q , 0.42), which were lower than the H e for all six 
units. This suggests non-random mating and/or popula- 
tion subdivision in each unit or the existence of null 
alleles. Among the six units, H e and H Q were the highest 
in the North-Central Alps (NCA) unit {H e = 0.85) and 
the Echigo-Mikuni (EM) unit (H Q = 0.69), respectively, 
while those of Western Chugoku (WC) were the lowest 
(H e = 0.33, H 0 = 0.19). In the NCA unit, the H Q value 
(H Q = 0.33) was significantly low (P < 0.001) compared 
with the highest H e value (H e = 0.85). The observed and 
expected heterozygosities of bears in western Japan 
(H e = 0.54, H Q = 0.31) were lower than that of bears in 
eastern Japan (H e = 0.73, H Q = 0.41). Nucleotide diversity 
(7T S ) based on synonymous substitutions within a unit 
ranged from 0.005 to 0.024 among the six units, and the 
overall value for black bears was 0.016. The jt s value of 
the WC unit (0.005) was the lowest among the six units. 

We calculated four statistics of genetic differentiation 
(K S t> Fst> Hedrick s G' ST and Jost s D) for pairs of man- 
agement units by using the nucleotide sequences of 
DQB variants (Table 5). Among the six management 
units, most 7<T ST and F ST values indicated statistically sig- 
nificant levels of genetic differentiation, including those 
between the WC and other units (mean K S t, 0.12; mean 



FsT) 0.17), and the Gassan-Asahi Iide (GAI) and other 
units (mean K ST , 0.05; mean F ST , 0.08). With respect to 
Hedricks G' ST and Jost s A these values trended to be 
high between the WC and other units. The high values 
of G st and D (more than 0.05) and significant values of 
K S t and F ST were observed only between the WC and 
EM/NCA/Hakusan-Okumino (HO) units (mean K S t> 
0.13; mean F ST , 0.20; mean G ST1 0.12; mean A 0.07). We 
calculated average values for each of four statistics be- 
tween management unit pairs in western and eastern 
Japan. The results indicated that the values for western 
Japan (K ST = 0.05, F ST = 0.10, G' ST = 0.05, D = 0.03) 
were higher than those for eastern Japan (K S t = 0.02, 
F ST = 0.05, G sx = 0.02, D = 0.01). 

Discussion 

Estimations of trans-species polymorphism in the DQB 
variants of the Japanese black bear 

In this study, we detected 44 DQB variants in 185 Japanese 
black bears. These 44 nucleotide variants translated to give 
31 allotypes and one putative, and the 31 allotypes were 
used to construct an ML tree to examine relationships 
among MHC class II sequences of ursine species (Figure 5). 
The tree showed that the DQB sequences of the Japanese 
black bear formed a monophyletic group with those of the 
giant panda. In our previous study, we reported that black 
bear DQB sequences cluster with DQB sequences of vari- 
ous mammalian species [39]. These results suggest that all 
MHC sequences of black bears are likely to have been 
derived from the DQB gene. 

The ML tree suggested a pattern of DRB trans-species 
evolution among ursine species, in agreement with the 



Table 4 Gene diversities and nucleotide diversities for each conservation and management unit 



MHC DQB exon 2 region (267 bp) mtDNA control region (Yasukochi et al. 2009) 





Conservation and 
management unit 


No. of 
samples 


No. of 
alleles 


H e 


H 0 




No. of 
samples 


No. of 
haplotypes 


h 


77 


Eastern 


Gassan-Asahi Iide (GAI) 


21 


4 


0.48 ± 0.08 


0.33 


0.018 ±0.009 


6 


5 


0.93 ± 0.12 


0.01 1 ± 0.008 


Japan 


Echigo-Mikuni (EM) 


13 


7 


0.79 ± 0.05 


0.69 


0.024 ±0.011 


16 


7 


0.82 ± 0.07 


0.009 ± 0.006 




North-Central Alps (NCA) 


12 


10 


0.85 ± 0.05* 


0.33* 


0.022 ±0.010 


23 


8 


0.85 ± 0.05 


0.011 ±0.007 d 




Hakusan-Okumino (HO) 


36 


18 


0.76± 0.05* 


0.39* 


0.019 ±0.008 


74 


8 


0.72 ± 0.03 


0.010 ± 0.006 


Western 


Northern Kinki (NK) 


23 


8 


0.68 ± 0.06 


0.56 


0.016 ±0.008 


96 (64) 


4 


0.75 ± 0.03 


0.011 ± 0.007 


Japan 


Eastern Chugoku (EC) 


8 


7 


0.83±0.07* 


0.38* 


0.013 ± 0.007 


31 (22) 


2 


0.25 ± 0.10 


0.003 ± 0.002 




Western Chugoku (WC) 


53 


11 


0.33 ± 0.06* 


0.19* 


0.005± 0.003 


154 (33) 


6 


0.56 ± 0.04 


0.005 ± 0.004 




Shikoku (SK) 


9 


6 


0.73 ± 0.10* 


0.33* 


0.021 ± 0.010 


11 c 


2 


0.55 ± 0.07 


0.009 ± 0.006 




Overall 


175 


44 


0.65 ± 0.03* 


0.36* 


0.016 ±0.007 


455 


24 


0.90 ± 0.01 


0.013 ± 0.008 



Bold character indicates conservation and management units which have more than 10 samples. Chokai Sanchi, Southern Ouu, Southern Alps and Kii Hanto units 
are excluded due to small number of samples. H e : the expected heterozygosity; H Q : the observed heterozygosity. 77 S or n: the nucleotide diversity, h: the haplotype 
diversity. *: a significant departure from Hardy-Weinberg equilibrium (P < 0.01). ( ): the number of haplotype data analyzed by [5]. a : the number includes 
subhaplotypes. b : the number was taken from the result of an identification of individuals using microsatellite analysis [90] because samples were obtained by hair 
trap. c : the number is samples that are newly added to samples used in our previous study [8]. d : the n value described in [8] {n = 0.001) was incorrect, and 
0.01 was a correct value. 
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Table 5 Four statistics of genetic differentiation for DQB variants between the pairs of management units 



Eastern Japan 



Western Japan 



Conservation and 


Gassan-Asahi 


Echigo-Mikuni 


North-Central 


Hakusan- 


Northern Kinki 


Eastern Western 


management unit 


lide (GAI) 


(EM) 


Alps (NCA) 


Okumino (HO) 


(NK) 


Chugoku (EC) Chugoku (WC) 


Eastern Japan 


GAI 


/Csi/fsT/Hedrick's 
















G' ST / Jost's D 














EM 


0.0407 0.0687 
















0.036**1 0.022** 














NCA 


0.057**/ 0.112**/ 


-0.006/ -0.01 2/ 














0.049**/ 0.030** 


-0.022*7-0.015** 












HO 


0.030**/ 0.063**/ 


0.008/ 0.029*/ 


-0.001/ 0.026*/ 












0.078*/ 0.049* 


-0.002*/ -0.002** 


-0.013*7 -0.010** 








Western Japan 


NK 


0.057**/ 0.104**/ 


0.017/0.054*/ 


-0.002/ 0.069*7 


-0.001/ 0.029*7 










0.063*/ 0.037* 


0.016**/ 0.011** 


0.021/0.015** 


0.023*7 0.016** 








EC 


0.009/ 0.116*7 


0.037*7 0.020/ 


0.029/ 0.031/ 


0.004/ 0.028*/ 


0.006/ 0.011/ 








0.0567 0.033** 


0.031*70.021** 


0.047*/ 0.033* 


0.062*/ 0.044* 


0.002*7 0.001** 






wc 


0.057*70.0737 


0.757*70.229*7 


0.141**/ 0.242**/ 


0.098**/ 0.118**/ 


0.140**/ 0.163**/ 


0.037*70.236*7 






0.026**/ 0.010** 


0.114/0.064 


0.114/0.065 


0.132/0.076 


0.105 / 0.055 


0.097/0.052 




SK 


0.000/ 0.032 / 


-0.010/ -0.008/ 


0.008/ -0.005/ 


0.005/ 0.008/ 


0.021/ 0.047*/ 


0.01 6/ 0.025/ 0.077**/ 0. 127**/ 






-0.003*7 -0.002** 


-o.oof7o.ooo** 


0.003*7 0.002** 


0.026*70.018** 


0.029*/ 0.018** 


0.025*7 0.01 6** 0.036**/ 0.017** 



From top line, the values in Table indicate the K ST and F ST , Hedrick's G' ST and Jost's D statistics. Bold italic shows values which are statistically significant for both 
of K ST and F STl with positive values of G' ST and D. Underline represents statistical significant K ST and F S t with more than 0.05 of G' ST and D. Bold characters indicate 
conservation and management units which have more than 10 samples. Chokai Sanchi, Southern Ouu, Southern Alps and Kii Hanto units are excluded due to 
small number of samples. *: P < 0.05; **: P < 0.01 (permutation test, 1,000 replicates). 



Goda et al [36] (Figure 5). The estimated divergence 
time of black bear DQB variants was approximately 11 ± 
4.6 MYA. Previous studies have estimated that a rapid 
radiation of the genus Ursus has occurred at around 5 to 
3 MYA [53,54]. Therefore, even if our estimate overesti- 
mated divergence time because of the limited number of 
synonymous sites, it is possible that the DQB allelic 
lineage of the Japanese black bear was shared with other 
species in the genus Ursus. In fact, Urth-DQBVIOI has 
been detected in the Asiatic black bear in Russia (YY, 
unpublished observation), although the PCR amplifica- 
tion indicating this was performed with the primer pair 
used in our previous study [39]. Fossil records suggest 
that the ancestral Japanese black bear migrated to the 
Japanese archipelago during the Middle Pleistocene (0.5 
MYA - 0.3 MYA) [55]. Therefore, in the Asiatic black 
bear, including the Japanese black bear, the Urth- 
DQB*01 lineage may have been maintained for at least 
0.3 million years (MY). 

If sequences of other ursine species are obtained in 
the future, it will be important to carefully distinguish 
orthologous genes from paralogous ones. It is possible 
that variants of paralogous genes can aberrantly reflect 
the trans-specific pattern of alleles in orthologous genes. 

Selective pressure on the DQB gene of the Japanese black 
bear 

To estimate whether balancing selection have acted 
on the DQB gene of black bears, we calculated d N /d s 
(or cj) ratios by using three methods: the modified 
Nei-Gojobori model, ML, and Bayesian inferences. 



All three methods indicated that the d N /d s values were 
significantly higher than unity, supporting the influence 
of balancing selection on the DQB locus of the Japanese 
black bear. According to the modified Nei-Gojobori 
method, the d N /d s value in the PBR was 3.8. This value 
is relatively higher than those of some other mamma- 
lian species (Table 2). The values of the giant panda 
DQB (1.3), the Namibian cheetah DRB (2.0), and the 
European mink DRB (1.8) are all lower than that of the 
black bear. In addition, the number of MHC variants in 
the black bear is considerably higher than that of the 
three species even though the variants of the Namibian 
cheetah may be derived from multiple loci (Table 2). 
Since these mammalian species are all designated as an 
endangered species in the International Union for Con- 
servation of Nature (IUCN) Red List [56], the effective 
population sizes of these species are expected to be 
quite small. In theory, the effect of selective pressure 
decreases with diminishing the effective population 
size, due to a larger effect of random genetic drift [57]. 
Thus, since the census population size of the Japanese 
black bear is relatively large, (more than 10,000 indivi- 
duals; the exact number of individuals is still unclear.), 
the black bear may show the higher d N /d s value, com- 
pared with other endangered species. 

Peptide-binding repertoires of the DQB molecules in the 
Japanese black bear 

The amino acids in the PBRs for the Japanese black bear 
DQBs can be divided into 15 different types (Figure 4). 
Theoretically, the average number of nonsynonymous 
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changes within the PBRs should be equal to the number 
of functionally different alleles that segregate in a popu- 
lation [58]. These groups might correspond to so-called 
"super-types" [59-61]. A super-type is a group of alleles 
with a similar repertoire of peptide binding. The mean 
number of nonsynonymous substitution in the PBRs 
(7<b) of the black bear DQB was 10.3, and this translated 
to an expected number of allelic lineages of approxi- 
mately ten. This is slightly inconsistent with the number 
of different types in the PBRs in Figure 4. This suggests 
that some PBRs may be functionally similar (i.e., a same 
super-type). Alternatively, since the length of DQB 
sequences is limited to exon 2 only, the resolution of the 
K B estimate may be lower. 

Effects of demographic events on the Japanese black 
bear population 

We compared results obtained from analyses of the 
mtDNA [8] and DQB genetic markers to estimate demo- 
graphic events that have occurred in Japanese black bear 
populations. MHC alleles are maintained for long peri- 
ods by balancing selection [62-64], while the coalescence 
times of mtDNA are short because of its maternal inher- 
itance mode and ploidy level. Due to this characteristic, 
mtDNA has a higher susceptibility to the effect of recent 
demographic factors. Therefore, the time scales in evolu- 
tionary history represented by these two types of genetic 
markers are different. 

The extent of sharing DQB allotypes among Japanese 
black bear management units was higher than that of 
mtDNA haplotypes (Figure 1). In particular, Urth- 
DQB*01 (Urth-DQB*0101) was frequently found 
throughout Japan. This observation could reflect the dif- 
ference in coalescence time between the two systems 
and a time-dependent difference in the degree of popu- 
lation fragmentation. Alternatively, if the same dominant 
pathogen affects different populations (i.e., a similar se- 
lective pressure), a DQB allotype that confers resistance 
to this pathogen is expected to be shared by different 
populations. On the other hand, the management unit- 
specific allotypes were found from most of the units. 
This suggests that rare allotypes also have been main- 
tained in the Japanese black bear population, possibly 
due to negative frequency-dependent selection [64]. 

Values of observed heterozygosity (H Q ) were lower 
than the H e for all six units. Some management units 
showed a significant departure from Hardy- Weinberg 
equilibrium (HWE) (Table 4). The excess of H e appeared 
not to be associated with the quality of DNA samples as 
there was no correlation between the number of homo- 
zygotes and the quality of DNA samples (P = 0.90). In 
addition, to examine the possibility of null alleles, we 
determined nucleotide sequences of the primer target 
region in 96 samples (YY, unpublished data). While all 



individuals showed only one degenerated site in the re- 
verse primer-annealing region, 10 individuals possessed 
seven degenerated sites in the forward primer region. In 
the forward primer region of the 10 individuals, there 
were two sets of nucleotide sequences that differed by 
less than 7 bp. The seven degenerated sites are derived 
from these two nucleotide sequences. However, the for- 
ward primer used in the present study (URDQBex2-F4) 
is designed to specifically anneal to only one particular 
sequence (see Methods) in the primer region. Therefore, 
the seven degenerated sites are unlikely to be the cause 
of null allele. This suggests that the excess of H e is not 
caused by technical problems. Since cross-locus amplifi- 
cation potentially decreases the possibility of the H e ex- 
cess (i.e., it increases the H Q values because of false 
heterozygotes), the H e excess suggests that our designed 
primer amplifies a single DQB locus. Although we can- 
not completely exclude the possibility of null alleles, the 
H e excess is more likely caused by demographic events 
such as non-random mating and/or population subdiv- 
ision in a local population. 

In our previous study [8], F ST values for mtDNA hap- 
lotypes among the management unit pairs in western 
Japan were larger than those among the pairs in eastern 
Japan. Since mtDNA haplotype distributions reflect re- 
cent demographic events, the large F ST in western Japan 
indicates a more infrequent gene flow within western 
Japanese populations. The level of genetic differentiation 
at DQB among the units in western Japan was also high. 
This increase may be due to non-random mating and 
genetic drift in western Japan in recent history. Non- 
random mating and genetic drift may be caused by frag- 
mentation and isolation of habitat. 

The genealogical structure of MHC alleles under 
strong symmetric balancing selection is similar to that of 
neutral genes, with the exception of a difference in time 
scales [65]. To estimate the effect of population subdiv- 
ision in the Japanese black bear in the distant past, we 
calculated the ratio of the largest genetic distance 
(TMRCA, ^tmrca) to the mean distance (d s or d) in 
bear DQB variants. At neutral genes, d TMRC A is expected 
to be twice the value of d [66]. The ratio of d TMRCA to 
d s in the DQB variants was approximately 3 (^tmrca = 
0.078 ± 0.032, d s = 0.028 ± 0.013), suggesting an effect 
of population fragmentation. Thus, the Japanese black 
bear population might have been isolated and fragmen- 
ted not only in the recent past but also in the relatively 
distant past. 

Genetic degradation of black bears in the WC unit 

For the WC unit in western Japan, values of genetic di- 
versity, for not only the MHC DQBs (H e = 0.33, n = 
0.005) but also mtDNA (h = 0.56, n = 0.005; [8]), were 
significantly lower than those for the other management 
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units (Table 4). In addition, both F S t and K S t values cal- 
culated between the WC unit and other units were sig- 
nificantly high, and also Hedrick s G ' ST and JostT) tended 
to be higher than those among the other units. This sug- 
gests that genetic drift and/or non-random mating also 
influence MHC diversity in the WC unit, and it implies 
that black bears in the WC unit might be gradually 
diminishing in their potential to adapt to environmental 
adversity, such as the presence of various pathogens. 
Thus, even if bears in the WC unit currently have resist- 
ance against infectious disease, we need to design a con- 
servation plan to maintain MHC diversity in the bears in 
this unit. The habitat area of the Eastern Chugoku (SK), 
Kii Hanto (KH), and Shikoku (SK) units in western 
Japan are smaller than that of the WC unit. Since previ- 
ous studies that examined genetic diversities of black 
bears by using mtDNA and microsatellite markers have 
suggested a lower level of genetic diversities in the popu- 
lations of western Japan [6,8], MHC variation in these 
management units might be lower than that in the WC 
unit. Therefore, it is essential to investigate MHC vari- 
ation in these units further by using an increased num- 
ber of samples. 

The effects of natural selection and genetic drift on 
patterns of MHC variation 

We observed significant H e excess and genetic differ- 
entiation among the management units due to popula- 
tion fragmentation. Ohnishi et al. [6] examined genetic 
diversity of Japanese black bears by using 10 microsat- 
ellite makers. In the previous study, the deviation of 
H Q from H e was not statistically significant in all 
examined local populations by using the HWE test. 
However, microsatellite markers that were not used in 
the study of Ohnishi et al. [6] showed such a deviation 
in the Eastern Chugoku (EC) and WC units in western 
Japan [4]. In addition, the genetic diversities within 
these units were similar to that of the isolated popula- 
tion of the brown bear [6]. Sutton et al. [67] proposed 
that negative frequency-dependent selection may be 
predominant in pre-bottlenecked populations and that 
greater loss of MHC diversity compared to neutral 
genetic diversity is accelerated because of uneven allele 
distributions. Such skewed MHC allele frequencies can 
occur under negative frequency-dependent selection or 
joint action of this selection [67,68]. Indeed, numerous 
MHC rare variants of the Japanese black bear were 
observed in the present study. Therefore, Japanese 
black bears might have lost some DQB rare variants 
due to a bottleneck event in their evolutionary history. 

Comparison of two criteria for MHC genotyping 

Due to the possibility of PCR artifacts, criteria for MHC 
genotyping on the basis of partial sequences (e.g., exon 2 



only) remain a matter of debate. For instance, previous 
studies have identified MHC variants according to more 
conservative criteria: If a rare sequence is detected in 
one individual alone (heterozygote) and the sequence 
differs by less than 3 bp from a known sequence of the 
same PCR product, the sequence is considered to be an 
artifact introduced by PCR error [69,70]. Although we 
attempted to carefully identify DQB variants, the meth- 
ods might not be definitive. Therefore, we compared the 
results of analyses between the two criteria. According 
to our criteria, we detected 44 variants from 185 bears. 
Application of more conservative criteria reduced the 
number to 23 variants and to 179 bear samples. Never- 
theless, results of analyses such as estimates of d N /d s 
ratio, genetic differentiation among units and deviation 
from the HWE did not largely change even though the 
number of variants was different between two criteria. 
This suggests that our conclusion is not affected by the 
differences between criteria. 

Future studies 

In this study, we detected a maximum of 44 variants in- 
cluding one putative (23 at a minimum). However, there 
is a possibility that some of these variants are pseudo- 
genized (i.e., contain frameshift or nonsense mutations 
in other exons or promoter regions), although we con- 
firmed that the nucleotide sequence of exon 2 in a near 
full-length cDNA sequence was identical to that of the 
Urth-DQB*01 [39]. Since a majority of our samples was 
collected from tooth, the size of amplified DNA frag- 
ments was limited to the 270 bp of exon 2. In addition, 
the sample size of our study was not sufficient to investi- 
gate MHC diversity in each management unit. To evalu- 
ate our results precisely, further analysis based on the 
sequences of the entire coding region of DQB should be 
performed by using more samples. 

We designed two pairs of primers to amplify the bear 
DQB exon 2. When one of these pairs (URDQBex2-Fl 
and -Rl) was used, three different sequences were 
detected from an individual by TA cloning, suggesting 
that at least two DQB loci were amplified. Wan et al. 
[31] have reported that the giant panda genome has 
three DQB loci, two functional and one non-functional. 
Kuduk et al. [37] have also confirmed expression of two 
DQB genes in the brown bear. The PCR results of our 
study suggest that the Japanese black bear may also pos- 
sess at least two DQB loci in the genome (See Methods). 
It will be desirable to determine the genomic sequence 
across the wider MHC region in order to assure the spe- 
cificity of each MHC locus in the bear genome. 

Conclusions 

In the present study, we detected 44 MHC variants from 
185 Japanese black bears. The phylogenetic analysis 
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suggests that their sequences are derived from the DQB 
loci. The value of d N /d s ratio indicates that the bear DQB 
variants have been maintained by balancing selection. 
The estimation of TMRCA among their variants suggests 
that the DQB lineages arose before the origin of genus 
Ursus. At least, the Urth-DQB*0101 variant has been 
shared between Japanese and continental black bears for 
0.3 MY (before the migration of Asiatic black bear from 
the Asian continent into Japanese archipelago). However, 
DQB genotyping for the Japanese black bear showed con- 
siderably lower H Q values than H e values in each local 
population. One may hypothesize that skewed allelic dis- 
tribution due to negative frequency-dependent selection 
and the bottleneck or population fragmentation in each 
management unit might accelerate loss of DQB rare var- 
iants. Therefore, the results of our study imply that MHC 
variation in Japanese black bears might have declined 
due to demographic events such as habitat fragmentation 
and isolation at the population or species level. Neverthe- 
less, DQB genetic diversity in Japanese black bears 
appears to be relatively high compared to some other 
endangered mammalian species. This result suggests 
that for the present the Japanese black bears appear to 
retain more potential resistance against pathogens than 
other endangered mammalian species. To prevent fur- 
ther decline of potential resistance against pathogens, a 
conservation policy for the Japanese black bear should 
be designed to maintain MHC rare variants in each of 
management units. To achieve this objective, it is im- 
portant to stop further fragmentation and isolation of 
bear habitats. 

Methods 

Collection of samples from Japanese black bears 

In this study, we collected 307 Japanese black bear sam- 
ples throughout the Japanese archipelago. The samples 
were collected from individuals that were hunted or cap- 
tured for pest control. Of these samples, 185 were used 
for MHC genotyping. These 185 samples were obtained 
from 12 conservation and management units designated 
by Yoneda [42] (Figure 1) between 1970 and 2006. Local 
populations of the black bears in Japan are divided into 
19 conservation and management units, taking into ac- 
count the prevention of bear movement, such as rivers, 
railways and so on. The division of local populations in 
our study followed the divisions of the bear management 
units. Five samples were obtained from the Chokai San- 
chi (CS) unit; 21 from the Gassan-Asahi Iide (GAI) unit; 
one from the Southern Ouu (SO) unit; 13 from the 
Echigo-Mikuni (EM) unit; one from the Southern Alps 
(SA) unit; 12 from the North-Central Alps (NCA) unit; 
36 from the Hakusan-Okumino (HO) unit; 23 from the 
Northern Kinki (NK) unit; 8 from the Eastern Chugoku 
(EC) unit; 53 from the Western Chugoku (WC) unit; 3 



from the Kii Hanto (KH) unit; and 9 from the Shikoku 
(SK) unit (Figure 1). We obtained the following types 
of samples: muscle and liver tissue, 42 samples; blood, 
10; hair, 12; oral mucous membrane, 4; and tooth or 
mandibula bone, 117. Based on the mtDNA analysis 
(see above), there were three genetically distinct groups 
of Japanese black bears. However, the sample sizes of 
the management units in the Shikoku/Kii-hanto popu- 
lations were quite small in this study (three samples in 
the KH unit and nine in the SK unit). Therefore, for 
the purpose of discussion of population structure, we 
combined the Shikoku/Kii-hanto population with the 
western Japanese population, as is the case with our 
previous study [8]. 

Experimental procedures for MHC genotyping 

DNA extraction was performed using the procedure 
described by our previous study [8]. To amplify the en- 
tire exon 2 of the DQB in the Japanese black bear, we 
initially used a single primer set URDQBex2-Fl and 
URDQBex2-Rl [39]. However, when this primer set was 
used, three variants per individual were detected from 
several individuals (i.e., at least two loci were amplified). 
When nucleotide sequences of the three variants were 
compared, two of three variants showed nearly identical 
sequences. Therefore, for amplification of a single DQB 
locus, a pair of internal primers for introns 1 and 2, 
URDQBex2-F4 (5/-GAGGCTGTGGGTTGGCCTGAGGC 
GCGGA-3/) and URDQBex2-R4 (5/-CCGCGGGGCGCG 
AACGGCCTGGCTCA-3/), was designed to specifically 
amplify the sequences in intron 1 region of their two similar 
variants. The PCR amplification, identification of heterozy- 
gous or homozygous individuals by TA cloning, and se- 
quencing were performed as described by Yasukochi et al. 
[39]. Murray et al. [71] showed that the probability that 
both of two variant sequences in a heterozygote are not 
detected by sequencing five clones (assuming that every 
variant is equally amplified at an equal frequency) is very 
low (P = 0.06). In our study, all final PCR products (188) 
were cloned, and 8 to 22 clones per individual were 
sequenced to decrease the possibility of false homozygotes 
caused by amplification bias. In addition, to confirm se- 
quencing reliability, we compared the result of direct- 
sequencing at the polymorphic sites in uncloned PCR 
product with cloned sequences (see Additional File 2: 
Figure S2). Since such a direct-sequencing are derived 
from the PCR product amplified from an original gDNA, 
it decreases the possibility of error incurred during clon- 
ing or PCR of the plasmid insert after cloning. If the fre- 
quencies of two variants detected from a heterozygote 
were extremely different and a double peak corresponding 
to the correct nucleotides was not observed after direct se- 
quencing, that individual was excluded from our analyses 
(see Additional File 2: Figure S2). In addition, we used high 
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fidelity polymerase during all PCRs (Ex Taq Hot Start Ver- 
sion TaKaRa Bio Inc.). 

Phylogenetic relationships and population differentiation 

The nucleotide sequences were aligned and translated 
into amino acids by using MEGA v5.0 [72]. The NJ and 
ML trees were constructed on the basis of Jones- 
Taylor-Thornton (JTT) model with gamma distribution 
[73]. The NJ tree construction and estimation of the 
best fit substitution model were performed by using the 
same software. The ML tree was estimated by using the 
PhyML 3.0 [74]. Bootstrap analysis was performed 
(1,000 replications in the NJ tree and 500 replications 
in the ML tree). A network of MHC variants was con- 
structed with the Neighbor-Net method [75] by using 
Splits-tree4 ver. 4.11.3 [76]. A Bayesian phylogenetic 
analysis was conducted with MrBayes v3.2.1 [77]. The 
analysis was implemented considering 11,500,000 gen- 
erations and tree sampling every 100 generations. The 
first 28,750 trees were discarded as burn-in. The genetic 
distance was estimated by using the JTT model. The 
Bayesian and ML trees were visualized with FigTree ver- 
sion 1.3.1 (http://tree.bio.ed.ac.uk/software/figtree/, [78]). 
The expected heterozygosity (H e ) for each management 
unit and Wrights F-statistics (F ST ) between the pairs of 
management units were calculated using Arlequin v3.11 
[79]. In addition, we calculated K S t [80] that, like F ST , is 
an index of genetic differentiation, using the DnaSP 5.10 
[81]. F ST is based on frequencies at polymorphic sites, 
while K S t uses the average number of differences be- 
tween sequences. The F ST and 7<T S t are possibly mislead- 
ing regarding the differentiation when heterozygosity is 
high. Thus, we also calculated Hedricks G' ST [82] and 
Josts D [83] using the SMOGD 1.2.5 [84]. 

Measure of genetic diversity and distance among 
DQB variants 

The extent of amino acid variability was evaluated using 
a Wu-Kabat plot [85]. The mean numbers of nonsynon- 
ymous substitutions per nonsynonymous site (d N ) and 
that of synonymous substitutions per synonymous site 
(d s ) were calculated for excluding identical sequences of 
DQB variants. We estimated nucleotide diversity (7T S ) 
among all the sequences in each management unit or 
across the entire population. The divergence time of the 
Japanese black bear DQB variants was estimated using 
synonymous substitutions alone. We assumed that the 
maximum genetic distance at synonymous sites roughly 
corresponded to the TMRCA of DQBs. These genetic 
distance parameters (ji s , d N , and d s ) were calculated 
using the modified Nei-Gojobori method [47] with 
Jukes-Cantor correction [86] by using MEGA v5.0. The 
transition/transversion bias in the calculation was 



assumed as R = 1.58 by the ML method. Standard errors 
were calculated from 1,000 bootstrapping replicates. 

Estimation of natural selection for DQB variants 

Ratios of d N /d s were calculated for each of the putative 
PBRs, non-PBRs, and for the overall exon 2 region. The 
extent of selective pressure at the amino acid level was 
measured by the ratio of nonsynonymous to synonym- 
ous rates (o = d N /d s ). We estimated o ratios with the 
maximum likelihood method by using CODEML in the 
PAML program version 4.2b [50]. Several models that 
differed in their hypotheses concerning the statistical 
distributions of the o ratio [48,49] were considered in 
this study. Nested models can be compared in pairs by 
using a likelihood-ratio test [87]. A Bayesian approach 
(Bayes empirical Bayes) in CODEML is usually used to 
predict codons that were potentially under positive se- 
lection or balancing selection in models M2a and M8 
[49]. Although it is difficult to distinguish between posi- 
tive selection and balancing selection, we use this Bayes- 
ian approach to estimate sites under balancing selection 
because MHC genes are thought to be affected by balan- 
cing selection rather than positive selection [65,88]. The 
ML method has the potential to overestimate o values 
because of the existence of recombination [70]. There- 
fore, we estimated the o value by using a Bayesian 
method with OmegaMap [51], which enabled us to infer 
balancing selection in the presence of recombination. 
The average number of nonsynonymous substitutions in 
the PBRs (Kb) was also estimated using method II, as 
described by Satta et al. [89]. 

Test for reciprocal recombination or gene conversion 

The GENECONV program was used with default settings. 
The global test for reciprocal recombination or gene con- 
version events was performed with 10,000 permutations of 
the sequence alignment to assess significance. Moreover, 
we examined the possibility of recombination by using the 
method described by Satta [44]. In the method, it is 
assumed that the number of substitutions in a particular 
region to the entire region is binomially distributed. We 
divided exon 2 into two regions, sheet and a helix, for 
the examination of recombination. The RDP4 software was 
also used for detecting a possible recombination [45]. 

Additional files 



Additional file 1: The NJ and Bayesian trees of ursine DQB and DRB 
genes. The NJ and Bayesian trees were constructed based on amino acid 
sequences of partial exon 2 region of ursine DQB or DRB genes. The 
model with gamma distribution is used for the NJ and Bayesian trees, 
respectively. Numbers in parentheses are Genbank accession numbers. 
Only bootstrap values over 50% and posterior probabilities over 80% are 
shown in the NJ and Bayesian trees, respectively. 
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Additional file 2: The comparison of nucleotides at the 
polymorphic sites between uncloned and cloned PCR products. 

(A) An example that an individual is identified as a heterozygote. The 
upper figure is the alignment of nucleotide sequences of uncloned and 
cloned PCR products. The lower figure is the result of direct-sequencing 
for uncloned PCR product. (B) An example that an individual is not 
identified as a heterozygote. Such an individual is excluded from analyses 
of the present study. 
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